

%Run programs
clc
clear all
close all
% 
%Initial steady-state
run Parameters
run BGP_growth

ss_c=C;
ss_y=Y;
ss_yT=YT;
ss_yNT=YNT;
ss_ha=Ha;
ss_hr=H_r;
ss_g=g(length(g));
ss_TT=Tnew;
ss_TNT=TNT;
ss_pi=Pi;
ss_piT=PiT;
ss_piNT=PiNT;
ss_p=P;
ss_pT=PT;
ss_pNT=PNT;
ss_PishareT=pi_shareT;
ss_PishareNT=pi_shareNT;
ss_w=w;
ss_vinnovator=vinnovator;
ss_jinnovator=Jinnovator;
ss_v=vadopter;
ss_j=Jadopter;
ss_vinnov=vinnov;
ss_chi = chi;
ss_YL=Y_L;
ss_eps=eps;
ss_Z=ZZ;
ss_AT=A_Z;
ss_rp=roy;
ss_IM=LHS;
ss_EX=RHS;
ss_IMT=IM_T;
ss_EXT=EX_T;
ss_IMNT=IM_NT;
ss_EXNT=EX_NT;
ss_Yw=Yw;
ss_Deltaeps=Delta_eps;
ss_rhochi = rho_chi; 

save Initial_values.mat ss_rhochi ss_chi ss_Yw Q ss_TNT ss_y ss_yT ss_yNT ss_pT ss_pNT ss_IMT ss_EXT ss_IMNT ss_EXNT ss_Deltaeps ss_c ss_y ss_ha ss_hr ss_g ss_TT ss_pi ss_piT ss_piNT ss_p ss_PishareT ss_PishareNT ss_w ss_vinnovator ss_jinnovator ss_v ss_j ss_vinnov ss_YL ss_eps ss_Z ss_AT ss_rp ss_EX ss_IM;
save Q.mat Q
save epsbar.mat epsbar
save lam.mat lam


%%Counterfactual steady-state
run Parameters_counterf
run BGP_growth_counterf


ss1_c=C;
ss1_y=Y;
ss1_yT=YT;
ss1_yNT=YNT;
ss1_ha=Ha;
ss1_hr=H_r;
ss1_g=g(length(g));
ss1_TT=Tnew;
ss1_TNT=TNT;
ss1_pi=Pi;
ss1_piT=PiT;
ss1_piNT=PiNT;
ss1_p=P;
ss1_pT=PT;
ss1_pNT=PNT;
ss1_PishareT=pi_shareT;
ss1_PishareNT=pi_shareNT;
ss1_w=w;
ss1_vinnovator=vinnovator;
ss1_chi = chi;
ss1_jinnovator=Jinnovator;
ss1_v=vadopter;
ss1_j=Jadopter;
ss1_vinnov=vinnov;
ss1_YL=Y_L;
ss1_eps=eps;
ss1_Z=ZZ;
ss1_AT=A_Z;
ss1_rp=roy;
ss1_IM=LHS;
ss1_EX=RHS;
ss1_IMT=IM_T;
ss1_EXT=EX_T;
ss1_IMNT=IM_NT;
ss1_EXNT=EX_NT;
ss1_Yw=Yw;
ss1_Deltaeps=Delta_eps;
ss1_rhochi = rho_chi; 

save Final_values.mat ss1_rhochi ss1_chi ss1_Yw ss1_TNT ss1_y ss1_yT ss1_yNT ss1_pT ss1_pNT ss1_IMT ss1_EXT ss1_IMNT ss1_EXNT ss1_Deltaeps ss1_c ss1_y ss1_ha ss1_hr ss1_g ss1_TT ss1_pi ss1_piT ss1_piNT ss1_p ss1_PishareT ss1_PishareNT ss1_w ss1_vinnovator ss1_jinnovator ss1_v ss1_j ss1_vinnov ss1_YL ss1_eps ss1_Z ss1_AT ss1_rp ss1_EX ss1_IM;
save epsbar.mat epsbar

 clear all

  dynare 'Dynamic_GFT_exponents_anticipated.mod' savemacro 
  save DynamicGFT.mat

%%Compute statistics

load Initial_values.mat
load Final_values.mat

dYL=log(ss1_YL./ss_YL)'.*100+(ss1_g-ss_g)*100/4
dC=log((ss1_c./ss1_p')./(ss_c./ss_p')).*100+(ss1_g-ss_g)*100/4./(1-bet)
g0=ss_g
g1=ss1_g
dhts=log(diag(ss1_PishareT)./diag(ss_PishareT)).*100
dimport_share=log((1-diag(ss1_PishareT))./(1-diag(ss_PishareT))).*100
dhr=log(ss1_hr./ss_hr)'.*100;
dha_total=log(sum(ss1_ha')./sum(ss_ha'))'.*100;
dhr_y=log((ss1_hr'./ss1_y)./(ss_hr'./ss_y)).*100;
dT=log((ss1_TT./ss1_TT(3))./(ss_TT./ss_TT(3))).*100;
dha_y=log((sum(ss1_ha')./ss1_y')./(sum(ss_ha')./ss_y'))'.*100;
Tss0=ss_TT./ss_TT(M);
Tss1=ss1_TT./ss1_TT(M);
%Dynamic gains from trade

load DynamicGFT.mat
% 
TT=150;
betg=bet; 
kk=[1:TT];
t=[1:TT];
TT=150;
betg=bet; 
kk=[1:TT];

for i=1:1

% Uold1(i)=sum((betg.^([1:TT]-1))'.*log(exp(resultsC1(1))))+sum((betg.^([1:TT]-1))'.*(cumsum(ones(TT,1).*(1+resultsgc1(1)))));
% Unew1(i)=sum((betg.^([1:TT]-1))'.*log(exp(resultsC1(1:TT)')))+sum((betg.^([1:TT]-1))'.*(cumsum(ones(TT,1).*((1+resultsgc1(1:TT)'./(1))))));
% 



Uold1=sum((betg.^([1:TT]).*log((1+gc1(1)).^[1:TT])));
prodgc1 = cumprod((1+gc1(1:TT)));
Unew1=sum((betg.^([1:TT]).*log(prodgc1)'));


%GFT1(i)=[(exp((1-betg)*(Unew1(i)-Uold1(i))))-1]*100;


% Uold3(i)=sum((betg.^([1:TT]-1))'.*log(exp(C3(1))))+sum((betg.^([1:TT]-1))'.*(cumsum(ones(TT,1).*(1+resultsgc3(1)))));
% Unew3(i)=sum((betg.^([1:TT]-1))'.*log(exp(resultsC3(1:TT)')))+sum((betg.^([1:TT]-1))'.*(cumsum(ones(TT,1).*((1+resultsgc3(1:TT)'./(1))))));


Uold3=sum((betg.^([1:TT]).*log((1+gc3(1)).^[1:TT])));
prodgc3 = cumprod((1+gc3(1:TT)));
Unew3=sum((betg.^([1:TT]).*log(prodgc3)'));




GFT3_anticipated=[exp((1-betg)*((log(exp(C3(TT))*(1+gc3(TT)))-log(exp(C3(1))*(1+gc3(1))))/(1-betg)+Unew3(i)-Uold3(i)))-1]*100;
GFT1_anticipated=[exp((1-betg)*((log(exp(C1(TT))*(1+gc1(TT)))-log(exp(C1(1))*(1+gc1(1))))/(1-betg)+Unew1(i)-Uold1(i)))-1]*100;


GFT3_BGP_anticipated =  [exp((1-betg)*((log(exp(C3(TT))*(1+gc3(TT)))-log(exp(C3(1))*(1+gc3(1))))/(1-betg)+sum(t.*beta.^t.*log(1+gc3(TT)))-sum(t.*beta.^t.*log(1+gc3(1)))))-1]*100
GFT1_BGP_anticipated =  [exp((1-betg)*((log(exp(C1(TT))*(1+gc1(TT)))-log(exp(C1(1))*(1+gc1(1))))/(1-betg)+sum(t.*beta.^t.*log(1+gc1(TT)))-sum(t.*beta.^t.*log(1+gc1(1)))))-1]*100



end


save Results_NN.mat






plot(1:M_.maximum_lag+TT, (oo_.endo_simul(strmatch('gc1',M_.endo_names,'exact'),1:M_.maximum_lag+TT)))
plot(1:M_.maximum_lag+TT, (oo_.endo_simul(strmatch('gc3',M_.endo_names,'exact'),1:M_.maximum_lag+TT)))
%%plot(1:M_.maximum_lag+TT, exp(oo_.endo_simul(strmatch('Hr1',M_.endo_names,'exact'),1:M_.maximum_lag+TT))./exp(oo_.endo_simul(strmatch('Y1',M_.endo_names,'exact'),1:M_.maximum_lag+TT)))


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%Plot consumption
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

C1counterf_anticipated(:)=log(exp(C1(TT))*(1+gc1(TT)))-log(exp(C1(1))*(1+gc1(1)))+log((prodgc1'))-log((1+gc1(1)).^[1:TT]);
C3counterf_anticipated(:)=log(exp(C3(TT))*(1+gc3(TT)))-log(exp(C3(1))*(1+gc3(1)))+log((prodgc3'))-log((1+gc3(1)).^[1:TT]);

%%%%United States
figure
hold on 
graph_baseline=plot([zeros(1,10) C1counterf_anticipated(1:50).*100])
set(graph_baseline,'LineWidth',2, 'LineStyle', '-','color', 'b');
graph_zeros=plot(zeros(60,1))
set(graph_zeros,'LineWidth',2, 'LineStyle', '--','color', 'k');
legend('Baseline','Initial BGP')
%legend('Baseline','Only IPR','Initial BGP')

title(['\fontsize{16} United States'])
ylabel(['\fontsize{14} log consumption (rel. initial BGP trend)'])
xticks([1 6 11 16 21 26 31 36 41 46 51 56 61])
xticklabels({'-10','-5','0','5','10','15','20','25','30','35','40','45','50','55', '60'})

hold off



%%%%China
figure
hold on 
graph_baseline=plot([zeros(1,10) C3counterf_anticipated(1:50).*100])
set(graph_baseline,'LineWidth',2, 'LineStyle', '-','color', 'b');
graph_zeros=plot(zeros(60,1))
set(graph_zeros,'LineWidth',2, 'LineStyle', '--','color', 'k');
legend('Baseline','Initial BGP')
%legend('Baseline','Only IPR','Initial BGP')

title(['\fontsize{16} China'])
ylabel(['\fontsize{14} log consumption (rel. initial BGP trend)'])
xticks([1 6 11 16 21 26 31 36 41 46 51 56 61])
xticklabels({'-10','-5','0','5','10','15','20','25','30','35','40','45','50','55', '60'})

hold off


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%Figures mechanism
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


T=50;
C1_counterf_anticipated=[(exp(C1(1:T)))];
C3_counterf_anticipated=(exp(C3(1:T)));
WP1_counterf_anticipated=((exp(Y_L1(1:T)).*exp(gt3(1:T)/(sig-1))));
WP3_counterf_anticipated=((exp(Y_L3(1:T)).*exp(gt3(1:T)/(sig-1)))); 
YL1_counterf_anticipated=((exp(Y1(1:T)).*exp(gt3(1:T)/(sig-1))));
YL3_counterf_anticipated=((exp(Y3(1:T)).*exp(gt3(1:T)/(sig-1))));
RP1paid_counterf_anticipated=exp(rp13(1:T))./exp(Y1(1:T));
RP1received_counterf_anticipated=exp(rp31(1:T))./exp(Y3(1:T));
EM31_counterf_anticipated=A31(1:T)./Z1(1:T);
EM13_counterf_anticipated=A13(1:T)./Z3(1:T);
Ha_Y1_counterf_anticipated=(exp(Ha11(1:T))+exp(Ha12(1:T))+exp(Ha13(1:T)))./exp(YT1(1:T));
Ha_Y3_counterf_anticipated=(exp(Ha31(1:T))+exp(Ha32(1:T))+exp(Ha33(1:T)))./exp(YT3(1:T));
Hr_Y1_counterf_anticipated=(exp(Hr1(1:T))./exp(YT1(1:T)));
Hr_Y3_counterf_anticipated=(exp(Hr3(1:T))./exp(YT3(1:T)));
AT11_counterf_anticipated=exp(A11(1:T))./exp(T1(1:T));
AT33_counterf_anticipated=exp(A33(1:T))./exp(T3(1:T));
T1T3_counterf_anticipated=exp(T1(1:T));
V1_counterf_anticipated=exp(Vinnov1).*exp(gt3/(sig-1));
V3_counterf_anticipated=exp(Vinnov3).*exp(gt3/(sig-1));
T1_total_counterf_anticipated=exp(T1(1:T))./(exp(T1(1:T))+exp(T2(1:T))+exp(T3(1:T)));
T3_total_counterf_anticipated=exp(T3(1:T))./(exp(T1(1:T))+exp(T2(1:T))+exp(T3(1:T)));
Z1_total_counterf_anticipated=exp(Z1(1:T))./(exp(Z1(1:T))+exp(Z2(1:T))+exp(Z3(1:T)));
Z3_total_counterf_anticipated=exp(Z3(1:T))./(exp(Z1(1:T))+exp(Z2(1:T))+exp(Z3(1:T)));

gc1_counterf_anticipated=gc1(1:T);
gc3_counterf_anticipated=gc3(1:T);
gy1_counterf_anticipated=gy1(1:T);
gy3_counterf_anticipated=gy3(1:T);
gp1_counterf_anticipated=gp1(1:T);
gp3_counterf_anticipated=gp3(1:T);

gt1_counterf_anticipated=gt1(1:T);
gt3_counterf_anticipated=gt3(1:T);
YT1_Y_counterf_anticipated=exp(YT1(1:T))./exp(Y1(1:T));
YT3_Y_counterf_anticipated=exp(YT3(1:T))./exp(Y3(1:T));
HTS_1_counterf_anticipated=exp(pi_share11(1:T));
HTS_3_counterf_anticipated=exp(pi_share33(1:T));
HTS_2_counterf_anticipated=exp(pi_share22(1:T));

HTS_all_3_counterf_anticipated=((exp(pi_share33(1:T)).*exp(YT3(1:T))+exp(pi_shareNT33(1:T)).*exp(YNT3(1:T)))./(exp(Y3(1:T))))
HTS_all_1_counterf_anticipated=((exp(pi_share11(1:T)).*exp(YT1(1:T))+exp(pi_shareNT11(1:T)).*exp(YNT1(1:T)))./(exp(Y1(1:T))))
T3_TTotal_counterf_anticipated=exp(T3(1:T))./(exp(T1(1:T))+exp(T2(1:T))+exp(T3(1:T)));
TOT13_counterf_anticipated=exp(P1(1:T))./exp(P3(1:T));

figure

subplot(3,2,1)
graphha1=plot([1:T],(Ha_Y1_counterf_anticipated(1:T)).*100)
set(graphha1,'LineWidth',2);
ylabel('% relat. BGP')
title('Adoption intensity US')

subplot(3,2,2)
graphha3=plot([1:T],(Ha_Y3_counterf_anticipated(1:T)).*100)
set(graphha3,'LineWidth',2);
title('Adoption intensity China')
ylabel('% relat. BGP')


subplot(3,2,3)
graphhr1=plot([1:T],(Hr_Y1_counterf_anticipated(1:T)).*100)
set(graphhr1,'LineWidth',2);
title('R&D intensity US')
ylabel('% relat. BGP')



subplot(3,2,4)
graphhr3=plot([1:T],(Hr_Y3_counterf_anticipated(1:T)).*100)
set(graphhr3,'LineWidth',2);
title('R&D intensity China')
ylabel('% relat. BGP')





figure
subplot(1,2,1)
hold on
graphgt1=plot([1:150],((gc1(1:150)).*100))
set(graphgt1,'LineWidth',2);
graphgt3=plot([1:150],((gc3(1:150)).*100))
set(graphgt3,'LineWidth',2,'LineStyle','--');
legend('USA', 'China')
title('Productivity growth')
ylabel('Percent')
hold off
subplot(1,2,2)
plot((exp(Y1(1:150)))./L1./(exp(Y3(1:150))./L3),'LineWidth',2)
title('Relative Productivity US to China')






figure 

subplot(3,2,1)
graphrpp=plot([1:T],(RP1paid_counterf_anticipated(1:T)).*100);
set(graphrpp,'LineWidth',2);
title('Royalties paid to China')
ylabel('Percent')


subplot(3,2,2)
graphrpr=plot([1:T],(RP1received_counterf_anticipated(1:T)).*100);
set(graphrpr,'LineWidth',2);
title('Royalties paid by China')
ylabel('Percent')



subplot(3,2,3)
graphrpp=plot([1:T],(HTS_1_counterf_anticipated(1:T)).*100);
set(graphrpp,'LineWidth',2);
title('Home trade share US')
ylabel('Percent')

subplot(3,2,4)
graphrpr=plot([1:T],(HTS_3_counterf_anticipated(1:T)).*100);
set(graphrpr,'LineWidth',2);
title('Home trade share China')
ylabel('Percent')

subplot(3,2,5)
graphrpr=plot([1:T],(TOT13_counterf_anticipated(1:T)).*1);
set(graphrpr,'LineWidth',2);
title('Relative prices US vs China')

% subplot(3,2,6)
% graphrpr=plot([1:T],(TOTw13_counterf_anticipated(1:T)).*1);
% set(graphrpr,'LineWidth',2);
% title('Relative productivity US vs China')


figure 

subplot(2,2,1)
graphha1=plot([1:T],(Ha_Y1_counterf_anticipated(1:T)).*100)
set(graphha1,'LineWidth',2);
ylabel('Percent')
title('Adoption intensity US')

subplot(2,2,2)
graphha3=plot([1:T],(Ha_Y3_counterf_anticipated(1:T)).*100)
set(graphha3,'LineWidth',2);
title('Adoption intensity China')
ylabel('Percent')


subplot(2,2,3)
graphhr1=plot([1:T],(Hr_Y1_counterf_anticipated(1:T)).*100)
set(graphhr1,'LineWidth',2);
title('R&D intensity US')
ylabel('Percent')



subplot(2,2,4)
graphhr3=plot([1:T],(Hr_Y3_counterf_anticipated(1:T)).*100)
set(graphhr3,'LineWidth',2);
title('R&D intensity China')
ylabel('Percent')




figure 

subplot(2,2,1)
graphrpp=plot([1:T],(RP1paid_counterf_anticipated(1:T)).*100);
set(graphrpp,'LineWidth',2);
title('Royalties paid by the US (% GDP)')
ylabel('Percent')


subplot(2,2,2)
graphrpr=plot([1:T],(RP1received_counterf_anticipated(1:T)).*100);
set(graphrpr,'LineWidth',2);
title('Royalties paid by China (% GDP)')
ylabel('Percent')



subplot(2,2,3)
graphrpp=plot([1:T],(HTS_1_counterf_anticipated(1:T)).*100);
set(graphrpp,'LineWidth',2);
title('Home trade share US')
ylabel('Percent')

subplot(2,2,4)
graphrpr=plot([1:T],(HTS_3_counterf_anticipated(1:T)).*100);
set(graphrpr,'LineWidth',2);
title('Home trade share China')
ylabel('Percent')

save Anticipated GFT1_anticipated GFT3_anticipated GFT1_BGP_anticipated GFT3_BGP_anticipated C1counterf_anticipated C3counterf_anticipated C1_counterf_anticipated C3_counterf_anticipated WP1_counterf_anticipated WP3_counterf_anticipated YL1_counterf_anticipated YL3_counterf_anticipated RP1paid_counterf_anticipated RP1received_counterf_anticipated EM31_counterf_anticipated EM13_counterf_anticipated Ha_Y1_counterf_anticipated Ha_Y3_counterf_anticipated Hr_Y1_counterf_anticipated Hr_Y3_counterf_anticipated AT11_counterf_anticipated AT33_counterf_anticipated T1T3_counterf_anticipated V1_counterf_anticipated V3_counterf_anticipated T1_total_counterf_anticipated T3_total_counterf_anticipated Z1_total_counterf_anticipated Z3_total_counterf_anticipated gc1_counterf_anticipated gc3_counterf_anticipated gy1_counterf_anticipated gy3_counterf_anticipated gp1_counterf_anticipated gp3_counterf_anticipated gt1_counterf_anticipated gt3_counterf_anticipated YT1_Y_counterf_anticipated YT3_Y_counterf_anticipated HTS_1_counterf_anticipated HTS_3_counterf_anticipated HTS_2_counterf_anticipated HTS_all_3_counterf_anticipated HTS_all_1_counterf_anticipated T3_TTotal_counterf_anticipated TOT13_counterf_anticipated;


%%%%%%%%%%%%%%Decomposition welfare%%%%%%%%%%%%%%%%%%%%%%%%



for t=2:TT
    gyn1(t) = (exp(Y1(t))/exp(Y1(t-1)))+(1+gt3(TT-1))^(1/(sig-1))-2
    gyn3(t) = (exp(Y3(t))/exp(Y3(t-1)))+(1+gt3(TT-1))^(1/(sig-1))-2
    ghr1(t) = (exp(Hr1(t))/exp(Hr1(t-1)))+(1+gt3(TT-1))^(1/(sig-1))-2
    gha11(t) = (exp(Ha11(t))/exp(Ha11(t-1)))+(1+gt3(TT-1))^(1/(sig-1))-2
    gha12(t) = (exp(Ha12(t))/exp(Ha12(t-1)))+(1+gt3(TT-1))^(1/(sig-1))-2
    gha13(t) = (exp(Ha13(t))/exp(Ha13(t-1)))+(1+gt3(TT-1))^(1/(sig-1))-2
    ghr3(t) = (exp(Hr3(t))/exp(Hr3(t-1)))+(1+gt3(TT-1))^(1/(sig-1))-2
    gha31(t) = (exp(Ha31(t))/exp(Ha31(t-1)))+(1+gt3(TT-1))^(1/(sig-1))-2
    gha32(t) = (exp(Ha32(t))/exp(Ha32(t-1)))+(1+gt3(TT-1))^(1/(sig-1))-2
    gha33(t) = (exp(Ha33(t))/exp(Ha33(t-1)))+(1+gt3(TT-1))^(1/(sig-1))-2 
end


    gyn1(1) = gc1(1);
    gyn3(1) = gc1(1);
    ghr1(1) = gc1(1);
    gha11(1) = gc1(1);
    gha12(1) = gc1(1);
    gha13(1) = gc1(1);
    ghr3(1) = gc1(1);
    gha31(1) = gc1(1);
    gha32(1) = gc1(1);
    gha33(1) = gc1(1);


TT=150;




prodgy1 = cumprod((1+gyn1(1:TT)));
prodgy3 = cumprod((1+gyn3(1:TT)));



prodghr1 = cumprod((1+ghr1(1:TT)));
prodghr3 = cumprod((1+ghr3(1:TT)));

prodgha11 = cumprod((1+gha11(1:TT)));
prodgha12 = cumprod((1+gha12(1:TT)));
prodgha13 = cumprod((1+gha13(1:TT)));

prodgha31 = cumprod((1+gha31(1:TT)));
prodgha32 = cumprod((1+gha32(1:TT)));
prodgha33 = cumprod((1+gha33(1:TT)));



C1counterf(:)=log(exp(C1(TT))*(1+gc1(TT)))-log(exp(C1(1))*(1+gc1(1)))+log((prodgc1(1:TT)'))-log((1+gc1(1)).^[1:TT]);
C3counterf(:)=log(exp(C3(TT))*(1+gc3(TT)))-log(exp(C3(1))*(1+gc3(1)))+log((prodgc3(1:TT)'))-log((1+gc3(1)).^[1:TT]);

Y1counterf(:)=log(exp(Y1(TT))*(1+gyn1(TT)))-log(exp(Y1(1))*(1+gc1(1)))+log((prodgy1'))-log((1+gc1(1)).^[1:TT])';
Y3counterf(:)=log(exp(Y3(TT))*(1+gyn3(TT)))-log(exp(Y3(1))*(1+gc3(1)))+log((prodgy3'))-log((1+gc3(1)).^[1:TT])';

Hr1counterf(:)=log(exp(Hr1(TT))*(1+ghr1(TT)))-log(exp(Hr1(1))*(1+gc1(1)))+log((prodghr1'))-log((1+gc1(1)).^[1:TT])';
Ha11counterf(:)=log(exp(Ha11(TT))*(1+gha11(TT)))-log(exp(Ha11(1))*(1+gc1(1)))+log((prodgha11'))-log((1+gc1(1)).^[1:TT])';
Ha12counterf(:)=log(exp(Ha12(TT))*(1+gha12(TT)))-log(exp(Ha12(1))*(1+gc1(1)))+log((prodgha12'))-log((1+gc1(1)).^[1:TT])';
Ha13counterf(:)=log(exp(Ha13(TT))*(1+gha13(TT)))-log(exp(Ha13(1))*(1+gc1(1)))+log((prodgha13'))-log((1+gc1(1)).^[1:TT])';


Hr3counterf(:)=log(exp(Hr3(TT))*(1+ghr3(TT)))-log(exp(Hr3(1))*(1+ghr3(1)))+log((prodghr3'))-log((1+ghr3(1))'.^[1:TT])';
Ha31counterf(:)=log(exp(Ha31(TT))*(1+gha31(TT)))-log(exp(Ha31(1))*(1+gha31(1)))+log((prodgha31'))-log((1+gha31(1)).^[1:TT])';
Ha32counterf(:)=log(exp(Ha32(TT))*(1+gha32(TT)))-log(exp(Ha32(1))*(1+gha32(1)))+log((prodgha32'))-log((1+gha32(1)).^[1:TT])';
Ha33counterf(:)=log(exp(Ha33(TT))*(1+gha33(TT)))-log(exp(Ha33(1))*(1+gha33(1)))+log((prodgha33'))-log((1+gha33(1)).^[1:TT])';

Ha1counterf = Ha11counterf+Ha12counterf+Ha13counterf;
Ha3counterf = Ha31counterf+Ha32counterf+Ha33counterf;




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%Welfare decomposition





%%%%United States
figure
hold on 
graph_baseline=plot([zeros(1,10)  C1counterf(1:35)])
set(graph_baseline,'LineWidth',2, 'LineStyle', '-','color', 'b');
graph_baseline=plot([zeros(1,10) Y1counterf(1:35)])
set(graph_baseline,'LineWidth',2, 'LineStyle', '--','color', 'b');

graph_baseline=plot([ zeros(1,10) Hr1counterf(1:35)])
set(graph_baseline,'LineWidth',2, 'LineStyle', '--','color', 'r');

graph_baseline=plot([zeros(1,10) Ha11counterf(1:35)])
set(graph_baseline,'LineWidth',2, 'LineStyle', '-.','color', 'k');

graph_baseline=plot([zeros(1,10) Ha12counterf(1:35)])
set(graph_baseline,'LineWidth',2, 'LineStyle', '--','color', 'k');

graph_baseline=plot([zeros(1,10) Ha13counterf(1:35)])
set(graph_baseline,'LineWidth',2, 'LineStyle', '-','color', 'k');


graph_zeros=plot(zeros(45,1))
set(graph_zeros,'LineWidth',2, 'LineStyle', ':','color', 'k');
legend('C','Y','Hr','Ha11','Ha12','Ha13','Initial BGP')
%legend('Baseline','Only IPR','Initial BGP')

title(['\fontsize{16} United States'])
ylabel(['\fontsize{14} log consumption (rel. initial BGP trend)'])
xticks([1 6 11 16 21 26 31 36 41 46 51 56 61])
xticklabels({'-10','-5','0','5','10','15','20','25','30','35'})

hold off



%%%%China
figure
hold on 
graph_baseline=plot([zeros(1,10)  C3counterf(1:35)])
set(graph_baseline,'LineWidth',2, 'LineStyle', '-','color', 'b');
graph_baseline=plot([zeros(1,10) Y3counterf(1:35)])
set(graph_baseline,'LineWidth',2, 'LineStyle', '--','color', 'b');

graph_baseline=plot([ zeros(1,10) Hr3counterf(1:35)])
set(graph_baseline,'LineWidth',2, 'LineStyle', '--','color', 'r');

graph_baseline=plot([zeros(1,10) Ha31counterf(1:35)])
set(graph_baseline,'LineWidth',2, 'LineStyle', '-.','color', 'k');

graph_baseline=plot([zeros(1,10) Ha32counterf(1:35)])
set(graph_baseline,'LineWidth',2, 'LineStyle', '--','color', 'k');

graph_baseline=plot([zeros(1,10) Ha33counterf(1:35)])
set(graph_baseline,'LineWidth',2, 'LineStyle', '-','color', 'k');


graph_zeros=plot(zeros(45,1))
set(graph_zeros,'LineWidth',2, 'LineStyle', ':','color', 'k');
legend('C','Y','Hr','Ha31','Ha32','Ha33','Initial BGP')
%legend('Baseline','Only IPR','Initial BGP')

title(['\fontsize{16} China'])
ylabel(['\fontsize{14} log consumption (rel. initial BGP trend)'])
xticks([1 6 11 16 21 26 31 36 41 46 51 56 61])
xticklabels({'-10','-5','0','5','10','15','20','25','30','35'})

hold off








